function d = spatial_difference_x(s,dh)
% second order differencing
% indexing: s(j,i)

d = 0*s;

% i = 1, one sided
d(:, 1) = ( ...
    -3*s(:,1) ...
    +4*s(:,2) ...
    -s(:,3) ...
    ) / 2 / dh;
% interior
d(:, 2:end-1) = ( ...
    s(:, (2:end-1)+1) ...
    -2*s(:, 2:end-1) ...
    +s(:, (2:end-1)-1) ...
    ) / dh^2;
% i = end, one sided
d(:, end) = -( ...
    -3*s(:,end) ...
    +4*s(:,end-1) ...
    -s(:,end-2) ...
    ) / 2 / dh;
